/*
 * File:   main.cpp
 * Author: thiago
 *
 * Created on 8 de Novembro de 2009, 22:26
 */

//#define DEBUG true
//#define DEBUG_VALORES true
//#define DEBUG_OPENGL false

#include <stdlib.h>
#include <malloc.h>
#include <iostream>
#include <math.h>
#include <Particula.h>
#include "RadixSort/radixsort.h"
#include "Scan/Scan.h"

//Includes GL
#include <GL/glew.h>
#include <GL/glut.h>
#include <GL/glext.h>

//Includes CUDA
#include <cuda_runtime.h>
#include <cutil.h>
#include <vector_types.h>
#include <cuda_gl_interop.h>
#include <cutil_gl_error.h>
#include <assert.h>
#include "cudpp/cudpp.h"

using namespace std;

extern "C" void launch_CalculaIDGridCell(unsigned int blocksPerGrid, unsigned int threadsPerBlock,
        float3* posicoes, unsigned int* values, unsigned int* keys, int numeroParticulas,
        int numColunas, float deltaX, float deltaY, unsigned int* count, float XminGrid, float YminGrid);

extern "C" void launch_FindNeighbors(unsigned int blocksPerGrid, unsigned int threadsPerBlock,
        float3* posicoes, Particula* particulas, Particula* particulas2, unsigned int* values, unsigned int* keys,
        unsigned int* count, unsigned int* offSet, int numeroColunasGridCell, int numeroLinhasGridCell, int numeroParticulas,
        float deltaX, float deltaY, float XminGrid, float YminGrid, float massaEspecificaInicial, float massa);

extern "C" void launch_ZeraCountEOffSet(unsigned int blocksPerGrid, unsigned int threadsPerBlock, int numeroCelulas,
        unsigned int* count, unsigned int* offSet);

const float massa = 1;
const float massaEspecificaInicial = 1000;
//Vetores para controle das particulas.
float3* posicoesParticulas;
Particula* particulas;
Particula* particulasDevice;
Particula* particulas2;
Particula* particulasDevice2;
unsigned int* keys;
unsigned int* keysDevice;
unsigned int* values;
unsigned int* valuesDevice;
unsigned int* count;
unsigned int* offSet;
unsigned int* countDevice;
unsigned int* offSetDevice;
int numeroParticulas;
int numeroColunasGridCell;
int numeroLinhasGridCell;
int numeroCelulasGrid;
float XminGrid = 0.0, XmaxGrid = 100.0, YminGrid = 0.0, YmaxGrid = 100.0, deltaXGrid, deltaYGrid;
float recuoXGrid = 2.0, recuoYGrid = 2.0;
//float* deltaXGridDevice;
//float* deltaYGridDevice;
int threadsPerBlock, blocksPerGrid;
//Utilizado no desenho das particulas.
int colunas, linhas;
float posicaoX, posicaoY, passoX, passoY;
bool mostraGrid, mostraEixo;
//VBO
GLuint positionsVBO;
float3* VBO;
GLuint positionsVBO2;
float3* VBO2;
bool VBO_Ativo;
const bool gpu = true;

void desenhaEixos() {

    glBegin(GL_LINES);
    glColor3f(1.0, 0.0, 0.0);
    glVertex3f(1.0, 0.002, 0.0);
    glVertex3f(0.0, 0.002, 0.0);
    glVertex3f(1.0, 0.001, 0.0);
    glVertex3f(0.0, 0.001, 0.0);

    glColor3f(0.0, 1.0, 0.0);
    glVertex3f(0.001, 1.0, 0.0);
    glVertex3f(0.001, 0.0, 0.0);
    glVertex3f(0.000, 1.0, 0.0);
    glVertex3f(0.000, 0.0, 0.0);

    //		glColor3f(0.0, 0.0, 1.0);
    //		glVertex3f(0.0, 0.0, 1.0);
    //		glVertex3f(0.0, 0.0, 0.0);

    glEnd();

}

void desenhaGrid() {

    float posicaoI = YmaxGrid, posicaoJ = XminGrid;

    glBegin(GL_LINES);

    glColor3f(1.0, 0.4, 0.0);
    for (int i = 0; i <= numeroLinhasGridCell; i++) {
        glVertex3f(XminGrid, posicaoI, 0.0);
        glVertex3f(XmaxGrid, posicaoI, 0.0);
        posicaoI = posicaoI - deltaYGrid;
    }
    for (int i = 0; i <= numeroColunasGridCell; i++) {
        glVertex3f(posicaoJ, YmaxGrid, 0.0);
        glVertex3f(posicaoJ, YminGrid, 0.0);
        posicaoJ = posicaoJ + deltaXGrid;
    }

    glEnd();
}

//Constroi o VBO
void buildVBO(void) {

#ifdef DEBUG_GL
    cout << "CONSTRUINDO VBO 1!" << endl;
    cout << "ANTES DO GenBuffer!" << endl;
    glGenBuffers(1, &positionsVBO);
    cout << "ANTES DO BindBuffer com o positionVBO!" << endl;
    glBindBuffer(GL_ARRAY_BUFFER, positionsVBO);
    cout << "ANTES DO BufferData!" << endl;
    unsigned int size = numeroParticulas * sizeof (float3);
    glBufferData(GL_ARRAY_BUFFER, size, posicoesParticulas, GL_DYNAMIC_DRAW);
    cout << "ANTES DO BindBuffer com o 0!" << endl;
    glBindBuffer(GL_ARRAY_BUFFER, 0);

    cout << "CONSTRUINDO VBO 2!" << endl;
    cout << "ANTES DO GenBuffer2!" << endl;
    glGenBuffers(1, &positionsVBO2);
    cout << "ANTES DO BindBuffer com o positionVBO2!" << endl;
    glBindBuffer(GL_ARRAY_BUFFER, positionsVBO2);
    cout << "ANTES DO BufferData2!" << endl;
    unsigned int size2 = numeroParticulas * sizeof (float3);
    glBufferData(GL_ARRAY_BUFFER, size2, posicoesParticulas, GL_DYNAMIC_DRAW);
    cout << "ANTES DO BindBuffer com o 0!" << endl;
    glBindBuffer(GL_ARRAY_BUFFER, 0);
#else
    glGenBuffers(1, &positionsVBO);
    glBindBuffer(GL_ARRAY_BUFFER, positionsVBO);
    unsigned int size = numeroParticulas * sizeof (float3);
    glBufferData(GL_ARRAY_BUFFER, size, posicoesParticulas, GL_DYNAMIC_DRAW);
    glBindBuffer(GL_ARRAY_BUFFER, 0);

    glGenBuffers(1, &positionsVBO2);
    glBindBuffer(GL_ARRAY_BUFFER, positionsVBO2);
    unsigned int size2 = numeroParticulas * sizeof (float3);
    glBufferData(GL_ARRAY_BUFFER, size2, posicoesParticulas, GL_DYNAMIC_DRAW);
    glBindBuffer(GL_ARRAY_BUFFER, 0);
#endif

    CUDA_SAFE_CALL(cudaGLRegisterBufferObject(positionsVBO));

    CUDA_SAFE_CALL(cudaGLRegisterBufferObject(positionsVBO2));

    CUT_CHECK_ERROR_GL(); //cout << "glGetError(): " << glGetError() << endl;
}

// Função callback chamada para fazer o desenho

void display(void) {

    // Limpa a janela de visualização com a cor de fundo especificada
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    if (mostraEixo) {
        desenhaEixos();
    }

    if (mostraGrid) {
        desenhaGrid();
    }

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

#ifdef DEBUG_GL
    if (VBO_Ativo == true) {
        cout << "HABILITOU O GL_VERTEX_ARRAY" << endl;
        glEnableClientState(GL_VERTEX_ARRAY);
        // Especifica que a cor corrente é vermelha
        //         R     G     B
        glColor3f(0.0f, 0.6f, 1.0f);
        cout << "ANTES DO BindBuffer com o positionVBO!" << endl;
        glBindBuffer(GL_ARRAY_BUFFER, positionsVBO);
        cout << "ANTES DO GLVERTEXPOINTER" << endl;
        glVertexPointer(3, GL_FLOAT, 0, 0);
        cout << "ANTES DO GLDRAWARRAYS" << endl;
        glDrawArrays(GL_POINTS, 0, numeroParticulas);
        cout << "ANTES DO BindBuffer com o 0!" << endl;
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        glDisableClientState(GL_VERTEX_ARRAY);
    } else {
        cout << "HABILITOU O GL_VERTEX_ARRAY" << endl;
        glEnableClientState(GL_VERTEX_ARRAY);
        // Especifica que a cor corrente é vermelha
        //         R     G     B
        glColor3f(0.0f, 0.6f, 1.0f);
        cout << "ANTES DO BindBuffer com o positionVBO!" << endl;
        glBindBuffer(GL_ARRAY_BUFFER, positionsVBO2);
        cout << "ANTES DO GLVERTEXPOINTER" << endl;
        glVertexPointer(3, GL_FLOAT, 0, 0);
        cout << "ANTES DO GLDRAWARRAYS" << endl;
        glDrawArrays(GL_POINTS, 0, numeroParticulas);
        cout << "ANTES DO BindBuffer com o 0!" << endl;
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        glDisableClientState(GL_VERTEX_ARRAY);
    }
#else
    if (VBO_Ativo == true) {
        glEnableClientState(GL_VERTEX_ARRAY);
        // Especifica que a cor corrente é vermelha
        //         R     G     B
        glColor3f(0.0f, 0.6f, 1.0f);
        glBindBuffer(GL_ARRAY_BUFFER, positionsVBO);
        glVertexPointer(3, GL_FLOAT, 0, 0);
        glDrawArrays(GL_POINTS, 0, numeroParticulas);
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        glDisableClientState(GL_VERTEX_ARRAY);
    } else {
        glEnableClientState(GL_VERTEX_ARRAY);
        // Especifica que a cor corrente é vermelha
        //         R     G     B
        glColor3f(0.0f, 0.6f, 1.0f);
        glBindBuffer(GL_ARRAY_BUFFER, positionsVBO2);
        glVertexPointer(3, GL_FLOAT, 0, 0);
        glDrawArrays(GL_POINTS, 0, numeroParticulas);
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        glDisableClientState(GL_VERTEX_ARRAY);
    }
#endif
    // Executa os comandos OpenGL
    glFlush();

    glutSwapBuffers();

}

// Inicializa parâmetros de rendering

void init(void) {
    // Define a cor de fundo da janela de visualização como preta
    glClearColor(0.0f, 0.0f, 0.0f, 1.0f);
    /*
        glClearDepth(1.0);
        glEnable(GL_LIGHT0);
        glEnable(GL_LIGHTING);
        glEnable(GL_DEPTH_TEST);
     */

}

// Função callback chamada quando o tamanho da janela é alterado

void AlteraTamanhoJanela(GLsizei w, GLsizei h) {
    // Evita a divisao por zero
    if (h == 0) h = 1;

    // Especifica as dimensões da Viewport
    glViewport(0, 0, w, h);

    // Inicializa o sistema de coordenadas
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();

    // Estabelece a janela de seleção (left, right, bottom, top)
    if (w <= h)
        gluOrtho2D(0.0f, 100.0f, 0.0f, 100.0f/**h/w*/);
    else
        gluOrtho2D(0.0f, 100.0f/**w/h*/, 0.0f, 100.0f);
}

void printDeviceProps() {
    cudaDeviceProp prop;
    int dev;

    CUDA_SAFE_CALL(cudaGetDevice(&dev));
    cout << "dev = " << dev << endl;
    CUDA_SAFE_CALL(cudaGetDeviceProperties(&prop, dev));

    printf("Name:                    %s\n", prop.name);
    printf("Global Memory:           %d\n", prop.totalGlobalMem);
    printf("Shared Memory per Block: %d\n", prop.sharedMemPerBlock);
    printf("Register per Block:      %d\n", prop.regsPerBlock);
    printf("Warp Size:               %d\n", prop.warpSize);
    printf("Memory Pitch:            %d\n", prop.memPitch);
    printf("Max Threads Per Block:   %d\n", prop.maxThreadsPerBlock);
    printf("Max Threads Dimension:   %d %d %d\n", prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
    printf("Max Grid Size:           %d %d %d\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
    printf("Constant Memory:         %d\n", prop.totalConstMem);
    printf("Version:                 %d.%d\n", prop.major, prop.minor);
    printf("Clock Rate:              %d\n", prop.clockRate);
    printf("Texture Alignment:       %d\n", prop.textureAlignment);
}

void printPosicoesParticulas() {
#ifdef DEBUG
    for (int j = 0; j < numeroParticulas; j++) {
        cout << "A posição (X) " << j << " do vetor A e igual a " << posicoesParticulas[j].x << endl;
    }
//    getchar();

    for (int j = 0; j < numeroParticulas; j++) {
        cout << "A posição (Y) " << j << " do vetor A e igual a " << posicoesParticulas[j].y << endl;
    }
//    getchar();
#endif
}

void printDadosParticulas() {
#ifdef DEBUG
    for (int i = 0; i < numeroParticulas; i++) {
        cout << "id: " << values[i] << "\t idCell: " << keys[i] <<
                "\t Massa Esp: " << particulas[values[i]].massaEspecifica <<
                "\t Pressao: " << particulas[values[i]].pressao <<
                "\t Visc: " << particulas[values[i]].viscosidade << endl;
    }
//    getchar();
#endif
}

void printDadosParticulas2() {
#ifdef DEBUG
    for (int i = 0; i < numeroParticulas; i++) {
        cout << "id: " << values[i] << "\t idCell: " << keys[i] <<
                "\t Massa Esp: " << particulas2[values[i]].massaEspecifica <<
                "\t Pressao: " << particulas2[values[i]].pressao <<
                "\t Visc: " << particulas2[values[i]].viscosidade << endl;
    }
//    getchar();
#endif
}

void printCells() {
    for (int i = 0; i < numeroParticulas; i++) {
        cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
    }
//    getchar();
}

void printCount() {
    for (int i = 0; i < numeroCelulasGrid; i++) {
        cout << "count[" << i << "] : " << count[i] << endl;
    }
//    getchar();
}

void radixSort() {
    int keybits = 32;

    // Creat the RadixSort object
    RadixSort radixsort(numeroParticulas, false);
    radixsort.sort(keysDevice, valuesDevice, numeroParticulas, keybits);

}

void scan() {
    Scan scanner(numeroParticulas);
    scanner.CriaOffset(countDevice, offSetDevice, numeroParticulas);
}

void CalculaKernelsOnGPU() {
    if (VBO_Ativo) {
        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) & VBO, positionsVBO));

        launch_ZeraCountEOffSet(blocksPerGrid, threadsPerBlock, numeroCelulasGrid, countDevice, offSetDevice);

        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO, valuesDevice,
                keysDevice, numeroParticulas, numeroColunasGridCell, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);

        scan();

        // Copy result from device memory to host memory
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        //cudaMemcpy(values, valuesDevice, numeroParticulas*sizeof(unsigned int), cudaMemcpyDeviceToHost);
        CUDA_SAFE_CALL(cudaMemcpy(count, countDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(offSet, offSetDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        printCells();

        printCount();

        cout << "\n Fazendo o radixSort" << endl;
#endif

        radixSort();
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(values, valuesDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        for (int i = 0; i < numeroParticulas; i++) {
            if (i < numeroCelulasGrid) {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << "\t offSet: " << offSet[i]
                        << "\t count[" << i << "] : " << count[i] << endl;
            } else {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
            }
        }
//        getchar();
#endif

        launch_FindNeighbors(blocksPerGrid, threadsPerBlock, VBO, particulasDevice, particulasDevice2, valuesDevice, keysDevice,
                countDevice, offSetDevice, numeroColunasGridCell, numeroLinhasGridCell, numeroParticulas,
                deltaXGrid, deltaYGrid, XminGrid, YminGrid, massaEspecificaInicial, massa);

        //CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(particulas2, particulasDevice2, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        cout << "Imprimindo particulas t (particulas)." << endl;
        printDadosParticulas();
        cout << "\n Imprimindo particulas t + 1 (particulas2)." << endl;
        printDadosParticulas2();
#endif
        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO));

        VBO_Ativo = false;

    } else {

        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) & VBO, positionsVBO));

        launch_ZeraCountEOffSet(blocksPerGrid, threadsPerBlock, numeroCelulasGrid, countDevice, offSetDevice);

        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO, valuesDevice,
                keysDevice, numeroParticulas, numeroColunasGridCell, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);

        scan();

        // Copy result from device memory to host memory
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        //cudaMemcpy(values, valuesDevice, numeroParticulas*sizeof(unsigned int), cudaMemcpyDeviceToHost);
        CUDA_SAFE_CALL(cudaMemcpy(count, countDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(offSet, offSetDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        printCells();

        printCount();

        cout << "\n Fazendo o radixSort" << endl;
#endif

        radixSort();
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(values, valuesDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        for (int i = 0; i < numeroParticulas; i++) {
            if (i < numeroCelulasGrid) {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << "\t offSet: " << offSet[i]
                        << "\t count[" << i << "] : " << count[i] << endl;
            } else {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
            }
        }
//        getchar();
#endif

        launch_FindNeighbors(blocksPerGrid, threadsPerBlock, VBO, particulasDevice2, particulasDevice, valuesDevice, keysDevice,
                countDevice, offSetDevice, numeroColunasGridCell, numeroLinhasGridCell, numeroParticulas,
                deltaXGrid, deltaYGrid, XminGrid, YminGrid, massaEspecificaInicial, massa);

        //CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        cout << "Imprimindo particulas t + 1 (particulas)." << endl;
        printDadosParticulas();
        cout << "\n Imprimindo particulas t (particulas2)." << endl;
        printDadosParticulas2();
#endif
        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO));

        VBO_Ativo = true;
    }
}


void calculaKernels() {
    if (VBO_Ativo) {
        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) & VBO, positionsVBO));

        launch_ZeraCountEOffSet(blocksPerGrid, threadsPerBlock, numeroCelulasGrid, countDevice, offSetDevice);

        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO, valuesDevice,
                keysDevice, numeroParticulas, numeroColunasGridCell, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);

        scan();

        // Copy result from device memory to host memory
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        //cudaMemcpy(values, valuesDevice, numeroParticulas*sizeof(unsigned int), cudaMemcpyDeviceToHost);
        CUDA_SAFE_CALL(cudaMemcpy(count, countDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(offSet, offSetDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        printCells();

        printCount();

        cout << "\n Fazendo o radixSort" << endl;
#endif

        radixSort();
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(values, valuesDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        for (int i = 0; i < numeroParticulas; i++) {
            if (i < numeroCelulasGrid) {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << "\t offSet: " << offSet[i]
                        << "\t count[" << i << "] : " << count[i] << endl;
            } else {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
            }
        }
//        getchar();
#endif

        //CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(particulas2, particulasDevice2, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));

        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO));

        VBO_Ativo = false;

    } else {

        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) & VBO, positionsVBO));

        launch_ZeraCountEOffSet(blocksPerGrid, threadsPerBlock, numeroCelulasGrid, countDevice, offSetDevice);

        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO, valuesDevice,
                keysDevice, numeroParticulas, numeroColunasGridCell, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);

        scan();

        // Copy result from device memory to host memory
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        //cudaMemcpy(values, valuesDevice, numeroParticulas*sizeof(unsigned int), cudaMemcpyDeviceToHost);
        CUDA_SAFE_CALL(cudaMemcpy(count, countDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(offSet, offSetDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        printCells();

        printCount();

        cout << "\n Fazendo o radixSort" << endl;
#endif

        radixSort();
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(values, valuesDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        for (int i = 0; i < numeroParticulas; i++) {
            if (i < numeroCelulasGrid) {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << "\t offSet: " << offSet[i]
                        << "\t count[" << i << "] : " << count[i] << endl;
            } else {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
            }
        }
//        getchar();
#endif

        //CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));

        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO));

        VBO_Ativo = true;
    }

    for (int indiceValue = 0; indiceValue < numeroParticulas; indiceValue++) {
        int jAtual = 0;
        int iAtual = 0;
        int jMin = 0;
        int jMax = numeroColunasGridCell;
        int iMin = 0;
        int iMax = numeroLinhasGridCell;
        int idCell = 0;
        int kInicial = 0;
        int kFinal = 0;
        float h = deltaYGrid;
        if (deltaXGrid < deltaYGrid) {
            h = deltaXGrid;
        }


        if (indiceValue < numeroParticulas) {
            cout << "\n" << endl;
            int posicaoParticula = values[indiceValue];
            cout << "indiceValue: " << posicaoParticula << endl;
            jAtual = (posicoesParticulas[posicaoParticula].x / deltaXGrid);
            cout << "jAtual: " << jAtual << endl;
            iAtual = (posicoesParticulas[posicaoParticula].y / deltaYGrid);
            cout << "iAtual: " << iAtual << endl;
            //int idCellAtual = (iAtual * numeroColunasGridCell) + jAtual;

            //VERIFICAO DAS EXTREMIDADES.
            if (jAtual == 0) {
                jMin = 0;
            } else {
                jMin = jAtual - 1;
            }

            if (jAtual == (numeroColunasGridCell - 1)) {
                jMax = jAtual;
            } else {
                jMax = jAtual + 1;
            }

            if (iAtual == 0) {
                iMin = 0;
            } else {
                iMin = iAtual - 1;
            }

            if (iAtual == (numeroLinhasGridCell - 1)) {
                iMax = iAtual;
            } else {
                iMax = iAtual + 1;
            }

            cout << "jMax: " << jMax << endl;
            cout << "jMin: " << jMin << endl;
            cout << "iMax: " << iMax << endl;
            cout << "iMin: " << iMin << endl;

            if (VBO_Ativo) {
                particulas2[posicaoParticula].massaEspecifica = 0;
            } else {
                particulas[posicaoParticula].massaEspecifica  = 0;
            }

            for (int i = iMin; i <= iMax; i++) {
                for (int j = jMin; j <= jMax; j++) {
                    idCell = (i * numeroColunasGridCell) + j;

                    kInicial = offSet[idCell];
                    kFinal = offSet[idCell] + count[idCell];

                    for (int k = kInicial; k < kFinal; k++) {
                        int indiceParticula = values[k];
                        cout << "Indice Part: " << indiceParticula << endl;
                        float squarePart1 = (posicoesParticulas[indiceValue].x - posicoesParticulas[indiceParticula].x);
                        float squarePart2 = squarePart1 * squarePart1;
                        float squarePart3 = (posicoesParticulas[indiceValue].y - posicoesParticulas[indiceParticula].y);
                        float squarePart4 = squarePart3 * squarePart3;
                        float squarePart5 = squarePart2 + squarePart4;
                        float distancia = sqrt(squarePart5);
                        if (VBO_Ativo) {
                            particulas2[posicaoParticula].massaEspecifica += particulas[indiceParticula].massaEspecifica * (315 / (64 * 3.14 * pow(h, 9))) * pow((h*h - distancia*distancia), 3);
                        } else {
                            particulas[posicaoParticula].massaEspecifica += particulas2[indiceParticula].massaEspecifica * (315 / (64 * 3.14 * pow(h, 9))) * pow((h*h - distancia*distancia), 3);
                        }

                    }

                }
            }
        }
    }
//    glGenBuffers(1, &positionsVBO);
//    glBindBuffer(GL_ARRAY_BUFFER, positionsVBO);
//    unsigned int size = numeroParticulas * sizeof (float3);
//    glBufferData(GL_ARRAY_BUFFER, size, posicoesParticulas, GL_DYNAMIC_DRAW);
//    glBindBuffer(GL_ARRAY_BUFFER, 0);
}

void teclado(unsigned char key, int x, int y) {

    switch (key) {

        case 27	:	exit(0);
			break;
	
	
        case 'g':
        case 'G':
            if (mostraGrid) {
                mostraGrid = false;
#ifdef DEBUG
                cout << "MostraGrid = false;" << endl;
#endif
            } else {
                mostraGrid = true;
#ifdef DEBUG
                cout << "MostraGrid = true;" << endl;
#endif

            }
            break;
        case 'e':
        case 'E':
            if (mostraEixo) {
                mostraEixo = false;
#ifdef DEBUG
                cout << "MostraEixo = false;" << endl;
#endif

            } else {
                mostraEixo = true;
#ifdef DEBUG
                cout << "MostraEixo = true;" << endl;
#endif

            }
            break;
        case 's':
        case 'S':
            if(gpu){
                CalculaKernelsOnGPU();
            } else {
                calculaKernels();
                if(VBO_Ativo){
                    printDadosParticulas2();
                } else {
                    printDadosParticulas();
                }
            }
            break;
    }
            
    glutPostRedisplay();

}

void iniciaDados() {

    mostraGrid = true;
    mostraEixo = false;

    //deltaXGrid = (float *) malloc(sizeof (float));
    //deltaYGrid = (float *) malloc(sizeof (float));

    deltaXGrid = (XmaxGrid - XminGrid) / numeroColunasGridCell;
    deltaYGrid = (YmaxGrid - YminGrid) / numeroLinhasGridCell;

    numeroCelulasGrid = numeroColunasGridCell * numeroLinhasGridCell;
    count = (unsigned int *) malloc(numeroCelulasGrid * sizeof (unsigned int));
    offSet = (unsigned int *) malloc(numeroCelulasGrid * sizeof (unsigned int));
    //Aloca o vetor de posições na memoria do host
    posicoesParticulas = (float3 *) malloc(numeroParticulas * sizeof (float3));
    particulas = (Particula *) malloc(numeroParticulas * sizeof (Particula));
    particulas2 = (Particula *) malloc(numeroParticulas * sizeof (Particula));
    keys = (unsigned int *) malloc(numeroParticulas * sizeof (unsigned int));
    values = (unsigned int *) malloc(numeroParticulas * sizeof (unsigned int));

    linhas = numeroParticulas / colunas;

    passoX = ((XmaxGrid - recuoXGrid) - (XminGrid + recuoXGrid)) / colunas;
    passoY = ((YmaxGrid - recuoYGrid) - (YminGrid + recuoYGrid)) / linhas;

    posicaoX = XminGrid + recuoXGrid;
    posicaoY = YmaxGrid - recuoYGrid;

    //Inicia os dados das particulas e gera os dados.
    for (int i = 0; i < numeroParticulas; i++) {
        particulas[i].massaEspecifica = massaEspecificaInicial;
        particulas[i].viscosidade = 1;
        particulas[i].pressao = 1;
        particulas[i].velocidade = 0;
        particulas2[i].massaEspecifica = massaEspecificaInicial;
        particulas2[i].viscosidade = 1;
        particulas2[i].pressao = 1;
        particulas2[i].velocidade = 0;
        //Gera os values.
        values[i] = i;
    }

    threadsPerBlock = 512;
    blocksPerGrid = (numeroParticulas + threadsPerBlock - 1) / threadsPerBlock;
}

void printDadosIniciais() {
#ifdef DEBUG
    cout << "\n-------------------- DADOS DA SIMULACAO --------------------" << endl;
    cout << "Numero de particulas: " << numeroParticulas << endl;
    cout << "Numero de linhas do Grid: " << numeroLinhasGridCell << endl;
    cout << "Numero de colunas do Grid: " << numeroColunasGridCell << endl;
    cout << "Numero de colunas para distribuicao das particulas: " << colunas << endl;
    cout << "DeltaX: " << deltaXGrid << endl;
    cout << "DeltaY: " << deltaYGrid << endl;
    cout << "------------------------------------------------------------" << endl;
//    getchar();
#endif
}

void iniciaParticulas() {

    //PREENCHENDO AS POSIÇÕES X.
    int contElementosLinha = 0;
    for (int j = 0; j < numeroParticulas; j++) {
        posicoesParticulas[j].x = posicaoX;

        contElementosLinha++;
        posicaoX = posicaoX + passoX;
        if (contElementosLinha >= colunas) {
            posicaoX = XminGrid + recuoXGrid;
            contElementosLinha = 0;
        }
    }

    //PREENCHENDO AS POSIÇÕES Y.
    int contElementosColuna = 0;
    for (int j = 0; j < numeroParticulas; j++) {
        posicoesParticulas[j].y = posicaoY;
        contElementosColuna++;
        if (contElementosColuna >= colunas) {
            posicaoY = posicaoY - passoY;
            contElementosColuna = 0;
        }
    }

    //PREENCHENDO AS POSIÇÕES Z.
    for (int j = 0; j < numeroParticulas; j++) {
            posicoesParticulas[j].z = 0.0;
    }
}

void iniciaCuda() {

    //Alocando as constantes da aplicação.
    /*CUDA_SAFE_CALL(cudaMalloc((void **) &massa_d, sizeof (float)));
    CUDA_SAFE_CALL(cudaMemcpyToSymbol("massa_d", &massa, sizeof (float), 0, cudaMemcpyHostToDevice));
    CUDA_SAFE_CALL(cudaGetSymbolAddress((void **)&massa_d, "massa_d"));*/
    //Aloca o vetor de particulas na memoria do device
    CUDA_SAFE_CALL(cudaMalloc((void **) & particulasDevice, numeroParticulas * sizeof (Particula)));
    CUDA_SAFE_CALL(cudaMalloc((void **) & particulasDevice2, numeroParticulas * sizeof (Particula)));
    CUDA_SAFE_CALL(cudaMalloc((void **) & keysDevice, numeroParticulas * sizeof (unsigned int)));
    CUDA_SAFE_CALL(cudaMalloc((void **) & valuesDevice, numeroParticulas * sizeof (unsigned int)));
    CUDA_SAFE_CALL(cudaMalloc((void **) & countDevice, numeroCelulasGrid * sizeof (unsigned int)));
    CUDA_SAFE_CALL(cudaMalloc((void **) & offSetDevice, numeroCelulasGrid * sizeof (unsigned int)));
    //cudaMalloc((void **) & deltaXGridDevice, sizeof (float));
    //cudaMalloc((void **) & deltaYGridDevice, sizeof (float));

    CUDA_SAFE_CALL(cudaMemcpy(valuesDevice, values, numeroParticulas * sizeof (unsigned int), cudaMemcpyHostToDevice));
    CUDA_SAFE_CALL(cudaMemcpy(particulasDevice, particulas, numeroParticulas * sizeof (Particula), cudaMemcpyHostToDevice));
    CUDA_SAFE_CALL(cudaMemcpy(particulasDevice2, particulas2, numeroParticulas * sizeof (Particula), cudaMemcpyHostToDevice));
    //cudaMemcpy(keysDevice, keys, numeroParticulas*sizeof(unsigned int), cudaMemcpyHostToDevice);
    CUDA_SAFE_CALL(cudaMemcpy(countDevice, count, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyHostToDevice));
    //cudaMemcpy(deltaXGridDevice, deltaXGrid, sizeof (float), cudaMemcpyHostToDevice);
    //cudaMemcpy(deltaYGridDevice, deltaYGrid, sizeof (float), cudaMemcpyHostToDevice);


}

/*
 *
 */
int main(int argc, char** argv) {

    //Verifica se todos os parametros realmente foram passados na chamada do programa.
    if (argc < 5) {
        cout << "Falta parametros!" << endl;
        exit(EXIT_SUCCESS);
    }

    numeroParticulas = atof(argv[1]);
    numeroColunasGridCell = atof(argv[2]);
    numeroLinhasGridCell = atof(argv[3]);
    colunas = atof(argv[4]);

    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
    glutInitWindowSize(800, 600);
    glutInitWindowPosition(10, 10);
    glutCreateWindow("Simulador de fluido 2D");
    GLenum err = glewInit();
    if (GLEW_OK != err) {
        //Problem: glewInit failed, something is seriously wrong.
        fprintf(stderr, "Error: %s\n", glewGetErrorString(err));
    }

    //printDeviceProps();
    cudaGLSetGLDevice(0);

    iniciaDados();
#ifdef DEBUG
    printDadosIniciais();
    printDadosParticulas();
#endif

    iniciaParticulas();

#ifdef DEBUG
    printPosicoesParticulas();
#endif

#ifdef DEBUG_GL
    cout << "ANTES DO VBO!" << endl;
#endif
    buildVBO();
    VBO_Ativo = true;
#ifdef DEBUG_GL
    cout << "DEPOIS DO VBO!" << endl;
#endif

    init();
    iniciaCuda();

    //OPENGL
    glutKeyboardFunc(teclado);
    glutDisplayFunc(display);
    glutReshapeFunc(AlteraTamanhoJanela);
    //glutIdleFunc(Desenha);

    glutMainLoop();

    return (EXIT_SUCCESS);
}
